/* this will use both (q, ftNL) and separate the (h) and (f) components and try to fit the simulation data using B1, B2 provided.
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mfnfitboth.h"

int main(int argc, char *argv[]) {
  printf("USAGE: mfn q ftNL name FEEDER\n");
  double calA=8.20092e-10;
  // get these values from supplement.pdf which shows that the correction factor from the initial conditions code matches the theory predictions!
  double qfnl;
  qfnl = atof(argv[1])*atof(argv[2]);
  calA = calA * 0.64 / (0.64+8.1e-9*pow(qfnl,2.0));
  double chisq, chisqmin=2000., dcmin=10., B1min=10., B2min=10., mumin=1.0;
  float fNL=0;
  // compute chisq in 3 dimensional space (deltac, B1=Bh, B2=Bf)
  double dc, B1, B2, mu;
  dc = 1.65;
  B1= 1.39;
  mu=1.0;
  for (dc=1.5; dc < 1.7; dc=dc+0.01) {
  for (B2=0.0; B2<3.01; B2=B2+0.01) {
  	chisq = mfn(calA, dc, atof(argv[1]), atof(argv[2]), argv[3], B1, B2, atof(argv[4]));
  	//printf("chisq: %g\n", chisq);
  	if (chisq<chisqmin) {
	    chisqmin = chisq;
	    dcmin = dc;
  	  B1min=B1;
	    B2min=B2;
	    mumin=mu;
	    printf("dc: %g, B1: %g, B2: %g, mu: %g, chisq: %g\n", dc, B1, B2, mu, chisq);
	  } else {
//	    printf("dc: %g, B1: %g, B2: %g, mu: %g, chisq: %g\n", dc, B1, B2, mu, chisq);
}	 
  }
  }
 printf("NDF: %d\n", NDFall);

  // write 
  FILE *fp;
  char fname[150], Opbase[200];
 // sprintf(Opbase, "/home1/02539/sza5154/feedersim/Theory/Output1024");
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024"); // this may need modification
  double M;
  int z;
// generate nu M correspondence curves
  for (z=0; z<3; z++) {
	sprintf(fname, "%s/nuMz%d.txt", Opbase, z);
	fp =fopen(fname, "w");
 	for (M=1e13; M<4e15; M=M*1.1) {
		fprintf(fp, "%g %g\n", M, nuc(1.6, M, AC, 0, 0, z));
	}
}

  // generate some theory curves
  for (z=0; z<3; z++) {
    //printf("M: %g, nuc: %g, nucD: %g \n", 1e14,nuc(deltac, 1e14, A, q, ftNL, z), nucD(deltac, 1e14, A, q, ftNL, z));
    sprintf(fname, "%s/%s/mfn/theoryz%d.dat", Opbase, argv[3], z);
    fp = fopen(fname, "w");
    for (M=1e13; M<4e15; M=M*1.05) {
      fprintf(fp, "%g %g\n", M,ratio(M,calA,dcmin,atof(argv[1]),atof(argv[2]),z,B1min,B2min, mumin));
    }
    fclose (fp);
  }
printf("here\n");
}


double mfn(double A, double deltac, double q, double ftNL, char* name, double Bh, double Bf, double mu) {
  char fname[150], string[500];
  FILE *fp;
  FILE *fpfNL;
  double M, rat, ratfNL, ratM3;
  int which;
  int ndf,z, ndfall=0;
  // open the relevant file for reading
  char Opbase[200];
  //sprintf(Opbase, "/home1/02539/sza5154/feedersim/Theory/Output1024");
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024"); 
  double chi2allz=0;
  double nuc1;
  int nz1=0;
  for (z=0; z<3; z++) {
    sprintf(fname, "%s/%s/mfn/mbinz%d.dat.ratio.avg", Opbase, name, z);
    fp = fopen(fname, "r");
    float col1, col2, col4, chi2=0.;
    ndf = 0;
    while (!feof(fp)) {
      if (fgets(string, 499, fp)) {
	sscanf(string, "%g %g %*g %g", &col1, &col2, &col4);
	M = col1*9.648e11; // in h^-1 M_sun
 	nuc1=nuc(deltac,M, A, q, ftNL,z);
	if (col1>0 && nuc1 < mu) {
	  if (z==1) nz1++;
	  //if (M>2e14) { //try only the highest ones!
	  rat = ratio(M,A,deltac,q,ftNL,z,Bh,Bf,1.0);
	  chi2+=pow((rat-col2)/col4,2);
	  ndf++; //}
	  //printf("r3: %g, r:4 %g, r5: %g\n", ratio3(M, A, deltac, q, ftNL, z), ratio4(M, A, deltac, q, ftNL,z), ratio5(M, A, deltac, q,ftNL, z));
	}
      }
    }
   // printf("nz1: %d, \n", nz1);
    //printf("%s, z:%d, chi2: %g, ndf: %d, chi2/ndf: %g\n", name, z, chi2, ndf, chi2/ndf);
    fclose(fp);
    chi2allz+=chi2;
    ndfall+=ndf;
  }
  NDFall=ndfall;
  //printf("ndf: %d, nz1: %d \t ", ndfall, nz1);
  return chi2allz/(ndfall-2.0);
}

double ratio(double M, double A, double deltac, double q, double ftNL, double z, double Bh, double Bf, double mu) {
  return mu*(1.0+ratio3h(M, A, deltac, q, ftNL, z)+ratio3f(M, A, deltac, q, ftNL, z)+ratio4h(M,A,deltac,q,ftNL,z, Bh)+ratio4f(M,A,deltac,q,ftNL,z,Bf)+ratio5h(M,A,deltac,q,ftNL,z,Bh)+ratio5f(M,A,deltac,q,ftNL,z,Bf));
}

double ratio3h(double M, double A, double deltac, double q, double ftNL, double z) {
  return M3Rh(M,A,q,ftNL)*He3(nuc(deltac, M, A, q, ftNL,z))/6.0-M3RDh(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/(6*nucD(deltac,M,A,q,ftNL,z)) ;
}

double ratio3f(double M, double A, double deltac, double q, double ftNL, double z) {
  return M3Rf(M,A,q,ftNL)*He3(nuc(deltac, M, A, q, ftNL,z))/6.0-M3RDf(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/(6.0*nucD(deltac,M,A,q,ftNL,z))+ M2R(M, A, q, ftNL)*He2(nuc(deltac, M, A, q, ftNL, z)/2.0)-M2RD(M,A,q,ftNL)*He1(nuc(deltac,M,A,q,ftNL,z))/(2.0*nucD(deltac,M,A,q,ftNL,z)) ; // add M2R here so as to not affect the rescaling due to B1 and B2 for others.
}

double ratio4h(double M, double A, double deltac, double q, double ftNL, double z, double Bh) {
  return  M4Rh(M,A,q,ftNL,Bh)*He4(nuc(deltac,M,A,q,ftNL,z))/24-M4RDh(M,A,q,ftNL,Bh)*He3(nuc(deltac,M,A,q,ftNL,z))/(24.0*nucD(deltac,M,A,q,ftNL,z))+pow(M3Rh(M,A,q,ftNL),2.0)*He6(nuc(deltac,M,A,q,ftNL,z))/72.0-M3Rh(M,A,q,ftNL)*M3RDh(M,A,q,ftNL)*He5(nuc(deltac,M,A,q,ftNL,z))/(36.0*nucD(deltac,M,A,q,ftNL,z));
}

double ratio4f(double M, double A, double deltac, double q, double ftNL, double z, double Bf) {
  return  M4Rf(M,A,q,ftNL,Bf)*He4(nuc(deltac,M,A,q,ftNL,z))/24.0-M4RDf(M,A,q,ftNL,Bf)*He3(nuc(deltac,M,A,q,ftNL,z))/(24.0*nucD(deltac,M,A,q,ftNL,z));
}

double ratio5h(double M, double A, double deltac, double q, double ftNL, double z, double Bh)
{
  return M5Rh(M, A, q, ftNL,Bh)*He5(nuc(deltac, M, A, q, ftNL, z))/120.0-M5RDh(M,A,q,ftNL,Bh)*He4(nuc(deltac,M,A,q,ftNL,z))/(120.0*nucD(deltac,M,A,q,ftNL,z))+M3Rh(M,A,q,ftNL)*M4Rh(M,A,q,ftNL,Bh)*He7(nuc(deltac,M,A,q,ftNL,z))/(24.0*6.0)+pow(M3Rh(M,A,q,ftNL)/6.0,3.0)*He9(nuc(deltac,M,A,q,ftNL,z))/6.0;
}

double ratio5f(double M, double A, double deltac, double q, double ftNL, double z, double Bf)
{
  return M5Rf(M, A, q, ftNL,Bf)*He5(nuc(deltac, M, A, q, ftNL, z))/120.0-M5RDf(M,A,q,ftNL,Bf)*He4(nuc(deltac,M,A,q,ftNL,z))/(120.0*nucD(deltac,M,A,q,ftNL,z));
}

double nuc(double deltac, double M, double A, double q, double ftNL, double z) {
  double gf;
  if (z==0){gf = 0.76001;} else if (z==0.5) {gf = 0.59455;} else if (z==1) {gf = 0.473345;} else if (z==2) {gf = 0.327534;} else { gf=0.; printf("no growth function value!\n"); }
  return deltac*0.76001/(sqrt(d2R(M,A,q,ftNL))*gf);
}

double nucD(double deltac, double M, double A, double q, double ftNL, double z) {
  double gf;
  if (z==0){gf = 0.76001;} else if (z==0.5) {gf = 0.59455;} else if (z==1) {    gf = 0.473345;} else if (z==2) {gf = 0.327534;} else {  gf=0.; printf("no growth function value!\n"); }
  return -deltac*0.76001*(0.5*(A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
   pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),0.5)/d2R(M,A,q,ftNL)/gf;
}

double M2R(double M, double A, double q, double ftNL) {
  return d2R(M,A,q,ftNL)/d2R(M,AC,0,0)-1;
}

double M3Rh(double M, double A, double q, double ftNL) {
  return d3R1(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),1.5);
}

double M3Rf(double M, double A, double q, double ftNL) {
  return d3R2(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),1.5);
}

double M3RDh(double M, double A, double q, double ftNL) {
  return (-9.*pow(A,2)*ftNL*pow(q,2)*I3R1(M)*
      (A*I2R1D(M) + 
        pow(A,2)*pow(ftNL,2)*pow(q,2)*
         I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*
       pow(q,2)*I2R2(M),2.5) + 
   (6*pow(A,2)*ftNL*pow(q,2)*
      I3R1D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*
       pow(q,2)*I2R2(M),1.5);
}

double M3RDf(double M, double A, double q, double ftNL) {
  return (-12.*pow(A,3)*pow(ftNL,3)*pow(q,3)*I3R2(M)*
      (A*I2R1D(M) + 
        pow(A,2)*pow(ftNL,2)*pow(q,2)*
         I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*
       pow(q,2)*I2R2(M),2.5) + 
   (8*pow(A,3)*pow(ftNL,3)*pow(q,3)*
      I3R2D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*
       pow(q,2)*I2R2(M),1.5);
}

double M4Rh(double M, double A, double q, double ftNL, double Bh) {
  return An(4)*pow(Bh*M3Rh(M,A,q,ftNL)/6.0,2.0);
}

double M4Rf(double M, double A, double q, double ftNL, double Bf) {
  return Bn(4)*pow(Bf*M3Rf(M,A,q,ftNL)/8.0, 4.0/3.0);
}

double M4RDh(double M, double A, double q, double ftNL, double Bh) {
  return An(4)*2*Bh*M3Rh(M,A,q,ftNL)*Bh*M3RDh(M,A,q,ftNL)/36.0;
}

double M4RDf(double M, double A, double q, double ftNL, double Bf) {
  return Bn(4)*pow(Bf*M3Rf(M,A,q,ftNL), 1.0/3.0)*Bf*M3RDf(M,A,q,ftNL)/12.0;
}

double M5Rh(double M, double A, double q, double ftNL, double Bh) {
  return An(5)*pow(Bh*M3Rh(M,A,q,ftNL)/6.0, 3.0);
}

double M5Rf(double M, double A, double q, double ftNL, double Bf) {
  return Bn(5)* pow(Bf*M3Rf(M, A, q, ftNL)/8.0, 5.0/3.0);
}

double M5RDf(double M, double A, double q, double ftNL, double Bf) {
 return Bn(5)*Bf*M3RDf(M,A,q,ftNL)*pow(Bf*M3Rf(M,A,q,ftNL)/8.0,2.0/3.0)/12.0;
}

double M5RDh(double M, double A, double q, double ftNL, double Bh) {
  return An(5)*3.0*pow(Bh*M3Rh(M,A,q,ftNL)/6.0, 2.0)*Bh*M3RDh(M,A,q,ftNL)/6.0;
}

double M2RD(double M, double A, double q, double ftNL) {
  return -(((A*I2R1(M) + pow(A,2.0)*pow(ftNL,2.0)*pow(q,2.0)*I2R2(M))*I2R1D(M))/
      (AC*pow(I2R1(M),2))) + (A*I2R1D(M) + 
      pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M))/(AC*I2R1(M));
}

double M3RD(double M, double A, double q, double ftNL) {
  return A*A*ftNL*q*(-(12.0*A*(0.75*I3R1(M)+A*ftNL*ftNL*q*I3R2(M))*(I2R1D(M)+A*ftNL*ftNL*I2R2D(M)))/pow(q*q+A*I2R1(M)+A*A*ftNL*ftNL*I2R2(M),2.5)+(6.0*I3R1D(M)+8.0*A*ftNL*ftNL*q*I3R2D(M))/pow(q*q+A*I2R1(M)+A*A*ftNL*ftNL*I2R2(M),1.5)); //check this with MM setup.nb numerical results!
}

double An(double n) {return tgamma(n-1.0)*pow(2.0, n-3.);}
double Bn(double n) {return tgamma(n)*pow(2.,n-1.);}

double d3R(double M, double A, double q, double ftNL) {
  return d3R1(M, A, q, ftNL) + d3R2(M, A, q, ftNL);
}

double d3R1(double M,double A,double q,double ftNL){return 6.0*q*q*ftNL*A*A*I3R1(M);}
double d3R2(double M,double A,double q,double ftNL){return 8.0*pow(q*ftNL*A,3.0)*I3R2(M);}

double I3R1(double M) {
  return pow(M, -4.31865)*exp(373.355+88145.1/pow(log(M),2)-9265.59/log(M));
}

double I3R2(double M) {
  return pow(M, -7.48587)*exp(660.953+175586./pow(log(M),2)-17909./log(M));
}

double d2R(double M, double A, double q, double ftNL) {
  return d2R1(M, A, q, ftNL) + d2R2(M, A, q, ftNL);
}

double d2R1(double M,double A,double q,double ftNL){return A*I2R1(M);}
double d2R2(double M,double A,double q,double ftNL){return pow(q*ftNL*A,2)*I2R2(M);}

double I2R1(double M) {
  return pow(M, -2.83160)*exp(246.051+58591./pow(log(M),2)-6140.93/log(M));
}

double I2R2(double M) {
  return pow(M, -3.13172)*exp(271.811+64250.2/pow(log(M),2)-6746.78/log(M));
}
double I2R1D(double M) {
  return (-2.8315989340868626*exp(246.05052999073197 + 58590.98533491765/pow(log(M),2) - 6140.926322745515/log(M)))/pow(M,3.8315989340868626) + 
   (exp(246.05052999073197 + 58590.98533491765/pow(log(M),2) - 6140.926322745515/log(M))*(-117181.9706698353/(M*pow(log(M),3)) + 6140.926322745515/(M*pow(log(M),2))))/pow(M,2.8315989340868626);
}

double I2R2D(double M) {
  return (-3.1317214228523023*exp(271.81104558816355 + 64250.21701980513/pow(log(M),2) - 6746.78180491083/log(M)))/pow(M,4.131721422852302) + 
   (exp(271.81104558816355 + 64250.21701980513/pow(log(M),2) - 6746.78180491083/log(M))*(-128500.43403961026/(M*pow(log(M),3)) + 6746.78180491083/(M*pow(log(M),2))))/pow(M,3.1317214228523023);
}

double I3R1D(double M) {
  return (-4.318649684650349*exp(373.3549059620789 + 88145.12900844424/pow(log(M),2) - 9265.593589986682/log(M)))/pow(M,5.318649684650349) + 
   (exp(373.3549059620789 + 88145.12900844424/pow(log(M),2) - 9265.593589986682/log(M))*(-176290.2580168885/(M*pow(log(M),3)) + 9265.593589986682/(M*pow(log(M),2))))/pow(M,4.318649684650349);
}

double I3R2D(double M) {
  return (-7.485871213528042*exp(660.9532909886348 + 175585.91029280794/pow(log(M),2) - 17908.95641931269/log(M)))/pow(M,8.485871213528043) + 
   (exp(660.9532909886348 + 175585.91029280794/pow(log(M),2) - 17908.95641931269/log(M))*(-351171.8205856159/(M*pow(log(M),3)) + 17908.95641931269/(M*pow(log(M),2))))/pow(M,7.485871213528042);
}

// the Hermite polynomials

double He1(double x) {return x;}
double He2(double x) {return x*x-1.0;}
double He3(double x) {return pow(x,3.0)-3.0*x;}
double He4(double x) {return pow(x,4.0)-6.0*x*x+3.0;}
double He5(double x) {return pow(x,5.0)-10.*pow(x,3.)+15.*x;}
double He6(double x) {return pow(x,6.0)-15.*pow(x,4.)+45.*x*x-15.;}
double He7(double x) {return pow(x,7.0)-21.*pow(x,5.)+105.*pow(x,3.)-105.*x;}
double He8(double x) {return pow(x,8.0)-28.*pow(x,6.)+210.*pow(x,4.)-420.*pow(x,2.)+105.;}
double He9(double x) {return pow(x,9.0)-36*pow(x,7.0)+378.0*pow(x,5.0)-1260.0*pow(x,3.0)+945.0*x;}
